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The discrete Boltzmann equation for both the ideal and a non-ideal fluid is extended by adding 
Langevin noise terms in order to incorporate the effects of thermal fluctuations. After casting the 
CNl . fluctuating discrete Boltzmann equation in a form appropriate to the Onsager-Machlup theory of 

» ! ■ linear fluctuations, the statistical properties of the noise are determined by invoking a fluctuation- 

' dissipation theorem at the kinetic level. By integrating the fluctuating discrete Boltzmann equation, 

a fluctuating lattice Boltzmann equation is obtained, which provides an efficient way to solve the 
equations of fluctuating hydrodynamics for ideal and non-ideal fluids. Application of the framework 
to a generic force-based non-ideal fluid model leads to ideal gas-type thermal noise. Simulation 
\ results indicate proper thermalization of all degrees of freedom. 
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I. INTRODUCTION 



o . 

O , Thermal fluctuations arc important in a wide variety of mcsoscopic flows in soft matter and biological physics. For 
■ example, thermal fluctuations are needed to produce diffusion in colloidal suspensions, to excite internal degrees of 
freedom in polymers and membranes suspended in a fluid, and to generate capillary waves on a fluid- fluid interface 
r/2 . Even in a simple fluid, thermal fluctuations lead to subtle non-linear effects like long time tails and the renormalization 
of transport coefficients 0. Such effects are most pronounced near critical points where fluctuations dominate Q. 
*^ Theoretically, thermally fluctuating mesoscopic flows are most conveniently dealt with within the framework of 

^~f ' fluctuating hydrodynamics [1,0. This approach, pioneered by Landau and Lifshitz for simple fluids, adds stochastic 
terms to the mechanical equations describing the flow, thereby enabling its statistical mechanical description. An 
CNJ • important ingredient in this formulation is the fluctuation dissipation theorem (FDT) relating the variance of the 
k*" | stochastic terms to the Onsagcr coefficients of the fluid. In Landau and Lifshitz's original formulation, the FDT relates 
the variance of the random stresses to the viscosities of the fluid. Thus an FDT-respecting fluctuating hydrodynamic 
' description allows one to calculate correlations in equilibrium as well as responses which deviate only weakly from 
. equilibrium. This has led to the widespread application of the method to complex fluid flow (I. 

Numerically, the solution of fluctuating hydrodynamic equations poses a serious challenge 0, 8( ■ Numerical solutions 
, of the Navier-Stokes equations, particularly in complicated geometries, are difficult even in the absence of thermal 
fluctuations. This has led to a whole range of innovative numerical methods for the Navier-Stokes equations, of which 
the lattice Boltzmann method (LBM) stands out due to its simplicity and easily parallelizablc nature. The basic 
version of the method is used to model isothermal, incompressible fluid flow, at both high and low Reynolds numbers 
. Extensions of the basic method allow for flows of complex and multicomponent fluids with microstructure [To| . 

The basic LBM ignores thermal fluctuations. Indeed, LB methods developed out of the need to eliminate thermal 
fluctuations inherent in their ancestor, the latticegas method, which needed excessive averaging over the noise to 
extract results pertinent to mean fluid flow [lT| - [T3j . Thus, the LBM as a preaveraged version of the lattice gas 
excluded all thermal fluctuations, making it inapplicable to the whole range of complex fluid flows mentioned above. 
Thermal fluctuations were partially reinstated in the LBE in the seminal work of Ladd |l4|, where, closely following the 
formulation of Landau and Lifshitz, fluctuations arc added only to the transport degrees of freedom of the LBE, i.e. the 
stresses. However, besides the conserved moments of mass and momentum and the non-conserved transport moments 
representing the stresses in the fluid, there exist in any LBE model also higher-order degrees of freedom, the so-called 
ghosts. The ghost sector, which is coupled to the transport sector at all but the largest length scales, acts as a sink for 
the thermal stress fluctuations, and thereby destroys the balance between fluctuation and dissipation. Consequently, 
the thermalization of the fluid remains incomplete. This is reflected in a similarly incomplete thermalization of other 
degrees of freedom that may be coupled to the fluid, for example colloidal particles. It was shown in [r| [r| how 
a fluctuating lattice Boltzmann equation (FLBE) can be constructed which offers a fully consistent discretization 
of the equations of fluctuating nonlinear hydrodynamics for an isothermal ideal fluid. This method reinstates the 
fluctuations in the ghost sector that were absent in previous work, thereby achieving a complete thermalization of the 
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fluid at all length scales, not only the largest ones. Notably, the general theory of LB-Langcvin equations was also 
derived in an early account by jl7| |. 

Recently, the Langevin approach of [lf| [13], originally devised for the ideal gas, has been generalized to describe 
thermal fluctuations in the LBE for a class of non-ideal fluids There, non-ideal fluid models based on a square- 
gradient (Ginzburg-Landau) free energy functionals were considered, implying that density fluctuations in such a 
fluid are spatially correlated with a correlation length which is proportional to the theoretical width of the diffuse 
liquid-vapor interface. This is in contrast to the ideal gas, where equilibrium correlations are generally absent for 
all degrees of freedom. Clearly, the non-trivial structure factor of a non-ideal fluid has to be faithfully represented 
by the correlations of the LB population densities. Using results of continuum kinetic theory, a general ansatz for 
these correlations has been proposed, noticing, however, that certain models might require modifications to account 
for implementation-specific details [l8|. This was, in particular, found to be the case for the model of Swift et al. fl9| . 
where, owing to the underlying modified equilibrium distribution, a spatially correlated form of thermal noise arises. 

In the present paper, we show how thermal noise can be incorporated into both ideal and non-ideal fluid versions of 
the discrete- velocity Boltzmann equation (DBE). The DBE is a precursor of the LBE that arises from the continuum 
Boltzmann equation by restricting velocity space to a finite number of velocities, while keeping position space and 
time continuous 0, [2l|. The fluctuating DBE (FDBE) is put forward here as a unifying starting point to construct 
fluctuating LB models and as a theoretical link between two successful frameworks of non-equilibrium physics: the 
LB method and the theory of linear regression of fluctuations due to Qnsager and Machlup [22, [23| • The utility of the 
FDBE is illustrated by re-deriving the fluctuating LBEs of the ideal gas [l5[ and the modified-equilibrium model [l8| . 
Moreover, based on the model of He, Shan and Doolen [24[ , we propose a FDBE for a generic force-based non-ideal 
fluid model and outline the requirements that have to be fulfilled in order to obtain a mathematically well-defined 
noise covariance. It is shown that the noise of this model is of the same form as in the ideal gas, in agreement with the 
physical notion that reversible interactions, such as those present in this type of non-ideal fluid, should not contribute 
to dissipation. This derivation constitutes the central result of the present work. 

Due to the presence of lattice corrections in the LBE expressions of the densities and fluxes and the generally 
non-linear form of the advection operator, the analysis of force-based non-ideal fluid models is technically difficult 
within the FLBE approach presented in The FDBE constitutes in this case a preferable starting point for the 
analysis of thermal noise in the LBE. A word of caution, however, is in order, as the DBE necessarily ignores also 
details of the spatial discretization scheme of the force, which can be an important ingredient in a LB model [25j . 
Only as long as the discretization scheme is conservative and does not give rise to any additional irreversible terms, 
can both LBE and associated DBE be regarded as equivalent concerning their behavior under thermal noise. The 
present framework should thus be seen as complementary to the FLBE approach of [l8j ). 

To set our work in context, it must be mentioned that the study of fluctuations in the continuous Boltzmann 
equation has a long history. The Boltzmann stosszahlansatz effectively removes fluctuations from the Boltzmann 
equation, giving a mean-field description of the fluid. However, fluctuations can be restored by promoting the Boltz- 
mann equation into a Langevin equation. This idea was first implemented by Kadomtsev [261 ]. and followed by 
Abrikosov and Khalatnikov [27| . Bixon and Zwanzig [28j gave a heuristic derivation of the fluctuation-dissipation 
theorem for the fluctuating Boltzmann equation. The work of Fox and Uhlenbeck [29|, [3(3] cast the theory into the 
framework developed by Onsager and Machlup for the linear regression of fluctuations, suitably generalised for the 
mixed parity of the Boltzmann equation under time. Fox and Uhlenbeck paid adequate attention to the conservation 
laws of mass, momentum and energy when constructing the fluctuation terms and deriving the FDT relation for the 
variances. From this point of view, our work can be looked upon as a consistent discretization in velocity space of 
the Boltzmann-Langevin equation developed by the above authors and an extension thereof to include non-ideal fluid 
interactions. When combined with a numerical scheme for temporally integra ting the discrete- velocity Boltzmann- 
Langevin equation, we obtain, among others, the FLBE previously derived in [l5l Il7l. [l8j. Subsequent to the work 
of Fox and Uhlenbeck, Kac and Logan [3l[ showed that the fluctuating Boltzmann equation can be derived from a 



master equation description of fluctuations in phase space, due to Siegert [32j. The phenomcnological assumption of 
an additive noise could then be rigorously justified. In a recent work, Dunweg et al. (l6l | have obtained the ideal gas 
FLBE of 15] from a coarse-graining of a lattice- gas model. In this sense, their work is the analogue for the FLBE 
of the work of Kac and Logan for the fluctuating Boltzmann equation. Their work, based on deterministic equations 
for probability densities, is thus complementary to the approach to the ideal gas presented below, based on stochastic 
equations for the fluctuating variables. 

The paper is organised as follows. In the next section, we summarise the Langevin theory of fluctuations for the 
Boltzmann equation and present key results of kinetic theory necessary for the extension of the theory to the non-ideal 
fluid. In section III we discuss the theory of finite-dimensional moment expansions in the discrete velocity-space, which 
facilitates the construction of the noise and the derivation of the FDT. The main new contribution of the present 
work is discussed in section IV, where we introduce the fluctuating DBE for the ideal and non-ideal fluid models and 
provide expressions for the noise covariances. In section V, we show how the FDBE may be integrated in time to 
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obtain the FLBE, providing an efficient numericai method for soiving the equations of fluctuating hydrodynamics 
for ideal or non-ideal fluids. In section VI, simulations of a particular implementation of the force-based FLBE are 
performed for a number of basic test cases. We conclude with a discussion of further extensions of the method and 
applications to problems of complex fluid flow. 

II. FLUCTUATIONS IN THE BOLTZMANN EQUATION 

The phenomcnological theory of fluctuations in the Boltzmann equation of a dilute gas, as we mentioned before, 
has been studied by several authors. The treatment which is most relevant for us is that of Fox and Uhlenbeck 
[30I ]. They show that the linearized Boltzmann equation with a Langevin noise term is of the form considered in the 
Onsager-Machlup linear theory of fluctuations [22I [23^ , as summarized in the next section. 

A. Linear theory of fluctuations 

The theory of linear rcgresssion of fluctuations, as originally proposed by Onsager and Machlup [22|, [23|, treats 
fluctuations of variables which are either even or odd under time-reversal symmetry. The Boltzmann equation has a 
mixed character, due to the presence of the reversible advective term and the irreversible collision term (see below). Fox 
and Uhlenbeck, therefore, generalise the Onsager-Machlup theory to such situations [2{|. They consider fluctuations 
of a set of variables a(i) = [ai(t), 02 (i), . . . , aj\r(t)] whose probability distribution at equilibrium is given by 

F e "(a) = |cxp(^-ia.G- 1 .a^ (1) 

where Z is a normalization constant and at is the conjugate transpose of a. The matrix of the equal-time correlations 
of the a's is fixed by the entropy matrix G , 

G y = («i(*K(*)>> ( 2 ) 
where * denotes complex conjugation. The a's are taken to obey linear Langevin equations of the form 

^a = -L-a + £, (3) 

where, L is a general time-evolution operator whose eigenvalues must have non-negative real parts. By convention, 
L is often written in terms of a symmetric and antisymmetric matrix, L = S + A. Typically, the symmetric matrix 
S describes relaxation, and hence, dissipation, while the antisymmetric matrix A contains the reversible part of the 
dynamics (e.g. free streaming and thermodynamic interactions). The random noise term is a zero-mean Gaussian 
process, and is therefore specified entirely by its second moment. It can be shown pi. |29| that the Langevin equation 
generates trajectories whose one-point probability distribution converges to cq. ([1]) if and only if the noise variance 
satisfies a fluctuation-dissipation relation, 

(6(^(t')) = (LG + GLt)..5( i -t'). (4) 

The theoretical framework presented in this section provides the basis of our treatment of the fluctuating DBE. 

B. Fluctuating Boltzmann equation 

Before analysing fluctuations in the discrete Boltzmann equation, it is worthwhile to briefly review the application 
of the Onsager-Machlup theory to the continuum Boltzmann equation for a dilute gas as done by Fox and Uhlenbeck 
[30| . since both treatments are closely connected. Fox and Uhlenbeck show that the linearized Boltzmann equation 
with noise can be brought into the form of Eq.([3]) above, thus enabling the use of the fluctuation-dissipation relation 
(|4]) to obtain the variance of the noise term. These results can be immediately compared to the results of our own 



1 In the original treatment of Fox and Uhlenbeck, the inverse of G is called 'entropy matrix' 
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derivation of the noise for the case of the DBE. While the following treatment deals with the Boltzmann equation for 
an ideal gas, non-ideal gas interactions can in principle be introduced by an explicit mean-field force. The contribution 
of the interaction force to the linearized Langevin dynamics can then be accounted for by an additional reversible 
contribution to the time-evolution operator of the Langevin equation ([3]). However, rather than dealing with this 
situation for the continuum Boltzmann equation, we postpone this to section IIV1 where the discrete Boltzmann 
equation for a non-ideal gas is discussed. This allows us in the following to focus on the essential aspects of the 
treatment of the continuum Boltzmann equation. 

The Boltzmann distribution function /(r, c, t) represents the number of particles at time t in a volume element 
drdc around the point (r, c) in the six-dimensional phase space of the particles. For particles of mass /z in equilibrium 
at temperature T with a density p Q , this distribution is independent of time and the coordinate r, and has a global- 
Maxwellian form, 

*>^»(sfcr) i/2 »"(-2iST c1 )' (5) 

with d being the spatial dimension. Note that, in order to remain close to the conventions used in the context of the 
LBE, we define the distribution function in terms of the mass density po rather than the number density and take 
the velocity c as the fundamental phase-space variable instead of the momentum. Close to equilibrium, then, it is 
possible to expand the distribution function around its global Maxwellian form, 

f(v,c,t) = /(c) [l + /i(r,c,i)] 

and argue that h(r, c, i) satisfies the linearized Boltzmann equation, 

dth + c ■ V/i = - / dV/(c)A(c, c') h(r, c', t) (6) 



The collision term is written using the Hilbert-Enskog kernel A(c, c'). Hilbert and Enskog have shown that the kernel 
is symmetric in c and c', isotropic, and has non- negative eigenvalues. Additionally, from the d+ 2 conservation laws 
of mass, momentum and energy, it follows that the kernel has d + 2 null eigenvalues. The eigenspectrum of the kernel 
cannot, in general, be obtained analytically. However, the eigenfunctions form a complete basis for expansion of 
functions in velocity space. In what follows, only the knowledge of the null-space and the positive-definiteness of the 
remaining eigenvalues are needed. Notice that the first two terms in eq. are odd under time reversal, while the 
last term is even. It is this mixed character of the equation which necessitates Fox and Uhlenbecks generalization of 
the Onsager-Machlup formulation. 

To transform the linear Boltzmann equation into a form consistent with Eq.([3]), define 



a(r,c,t) = \J /(c) h(r, c, t) 



A(r,c;r',c') = J /(c) c ■ V<5(r - r') S(c - c') 



S(r,c;r',c') = ^/>)/(c')A(c, c')5(r - r') 

Treating the labels r and c as indices, A(r, c; r', c') is antisymmetric because of the factor V<5(r— r'), while S(r, c; r', c') 
is symmetric with non-negative eigenvalues because of the factor A(c, c'). In these variables, the linearized Boltzmann 
equation can be written as, 



d t a(r,c,t) + J d d r'd d c'A(r,c; r', c') a(r', c', t) 
= - J d d r'd d c'S{r, c; r', c') a(r', c', t) 



Noting the symmetry of the kernels, the above equation has the form of the standard regression equation of a Gaussian 
Markov process, 

|(a>+A-(a>=-S.(a) 

where A is antisymmetric in the indices, and S is symmetric in the indices and has non-negative eigenvalues. The 
phcnomenological theory of fluctuations in the Boltzmann equation is now reduced to adding noise terms which respect 
the conservation laws and reproduce Maxwell-Boltzmann equilibrium. 
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The linear regression equation of a(r, c, t) is transformed into a stochastic differential equation for a Gaussian 
Markov process by adding a noise term \J /(c) £(r, c, t). In the original variables, this equation gives, 

d t h + c • V/i = - J d d c'/(c)A(c, c') h(r, c', t) + £(r, c, t) 

where £(r, c, t) is a zero mean Gaussian random variable with a correlation that has to be determined from the FDT, 
cq. (UJ). The essential step is to expand the Boltzmann entropy to quadratic order in h(r, c, t) and obtain the entropy 
matrix G(r, c; r', c') = (a(r, c)a(r', c')) which occurs in the fluctuation-dissipation relation. Fox and Uhlenbeck obtain 

G- 1 (r, c; r', c') = S(r - r') *(c - c') , (7) 

reflecting the fact that, in a dilute gas, particles are essentially uncorrelatcd and obey Poissonian fluctuation statistics 
[lH . It follows immediately that 

(£(r, c, t) £(r', c', t')) = 2A(c, c') S(v - r') S(t - t') . (8) 

This is the fluctuation-dissipation relation for the linear fluctuating Boltzmann equation. It is intuitively clear why 
the variance of the fluctuation is determined by the collision kernel alone. The advection kernel contains the reversible 
part of the dynamics and only shifts the distribution function in phase space, without a change of shape. The collision 
kernel, on the other hand, causes dissipation and therefore a broadening of the distribution function. Thus, eq. flSJ) 
means that the strength of the fluctuations is related to the amount of dissipation - the essential content of the 
fluctuation-dissipation relation. 

An explicit construction of the noise is now possible in terms of the eigcnfunctions and eigenvalues of the collision 
operator. Assuming that there arc discrete eigcnfunctions <pi(c) with eigenvalues A^, which are orthonormal with 
respect to the weight function w(c) = f(c)/p 0l 

J dc 10(0)^(0)^(0) = Sij 

the kernel has the expansion 

00 

A(c J c') = 5^Ai&(c)&(c / ) 
»=i 

The noise can likewise be expanded in the basis of the eigenfunctions as 

00 

f(r,c,t) =Y^£i(r,t)4> i (c) 

i=l 

The conservation laws imply the collision kernel has d + 2 zero eigenvalues. Taking these (in three dimensions) to be 
i = 1, . . . , 5 we see that the summations must start from i = 6 in the expansion of the Hilbert-Enskog kernel, and in 
the construction of the noise, we must set £i(r, t) = for i = 1, . . . , 5, since the conserved modes do not contribute 
to dissipation. Then, it only remains to obtain the variances of the spatial coefficients of the noise and a little more 
algebra shows that these are determined by the eigenvalues of the collision kernel. The final result is 

(£i(r, t)^(r', f)) = 2 Aj 6(v - r') 5(t - t') . 

This result shows that every non-conserved eigenmode contributes to the dissipation, and therefore, every non- 
conserved mode must receive fluctuations from the noise to maintain Boltzmann equilibrium. The theory above 
provides a rigorous justification of the assumptions used intuitively in [l5j . 



C. Correlations in kinetic theory 

As seen from the above derivation, a crucial ingredient to the FDT of the Boltzmann equation is the equilibrium 
correlation matrix (entropy matrix) G = (h(r, c)h(r',c')) of the fluctuations h = (/ — /)// of the distribution function. 
This matrix can be determined if an expression for the equilibrium entropy in terms of the /i's is known, as is the 
case for the continuum Boltzmann equation due to the H-theorem. A more direct approach to G, that can moreover 
easily deal with non-ideal fluids, is based on the kinetic theory of fluctuations due to Klimontovich [l8|, HH HH • This 
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approach will form the basis for our treatment for fluctuations in the discrete Boltzmann equation for the non-ideal 
gas. 

The starting point is the notion of a phase-space density F, defined as 

F(r, c, t) = fi S ( r ~ r *W)% - Ci(t)) , 

i 

where r.; and denote the (time-dependent) positions and momenta of the individual particles of mass p of 
the fluid. As J F(r, c, t)drdc/fi = N, the quantity F(r, c, i)dvdc/p represents the instantaneous number of 
particles in the phase-space cell drdc. The average number of particles (times mass) (F) in a phase-space 
cell can be obtained by integrating F weighted by the full iV-particle distribution function f N , (F(r,c,t)) = 
J dri ■ ■ ■ dv^dci ■ ■ ■ dcjv/jv(ri, . . . , rjy, Ci, . . . cn)F(t, c, t). Obviously, this definition is equal to the one-particle dis- 
tribution function f\, 

(F(r,c,i)> =/i(r,c,t). 

The average of the second moment (F(r, c, t)F(r', c', t')) can be related to the two-particle distribution function 
/2(r, c, r', c', t, t') by separating from the double-sum of delta functions the part where the two particles are identical. 
This gives 

(F(r, c, t)F{v',c', t')) = / 2 (r, c, t, r',c', t') + pS(r - v')S{c - c')S(t - i')/i(r, c, f) . (9) 

Similar relations can be derived for all reduced n-particle distribution functions. In the present work, we will only 
need the reduced one- and two-particle distribution functions. 

The quantity /i(r, c, t) = (F(r, c, t)) denotes the ensemble-averaged number of particles (times mass) at r, c at time 
t. and consequently, its time evolution will approximately be governed by the deterministic Boltzmann equation [35j . 
Intuitively, this notion motivates us to interpret the quantity F itself as the instantaneous fluctuating one-particle 
distribution function // (r,c,t). However, as the fluctuations in F are due to collisions between a large number of 
particles, the assumption of molecular chaos suggests the time evolution of // to be approximately governed by the 
Boltzmann-Langevin equation 28 1. In this case, the fluctuations of // have Gaussian character with only their second 



moment being non-trivial, in contrast to the phase-space density F, which encapsulates also higher-order correlation 
functions. The essential relation for the correlations of fluctuations of the one-particle distribution function is provided 
by the ansatz 

(//(r,c,i)AV,c',0) = (F(r,c,i)F(r',c',0>, 

which, as we show below, allows us to make a connection via cq. ((9j) to a statistical mechanical model for the equal-time 
pair correlation function. 

For the purpose of deriving the FDT, one usually considers fluctuations around a global equilibrium state with 
density po and zero flow velocity, i.e. f\ in the previous definition is taken to be a global Maxwellian distribution, 
/i(r,c,i) = /i(c), with / given by eq. ([5]). The quantity <J/i(r, c,t) = // (r, c,t) — /i(c) then represents the fluctua- 
tions of the one-particle distribution function f\ over a uniform reference state described by /i(c). Restricting to a 
translationally invariant system, the equal-time correlation function follows as 



(S Mr, 0)6^(^,0')) = (F(r,c)F(r',c')> - /i(c)A(c') 

- /i(c)/i(c')[ff(r - r') - 1] + f ,5(r - v')S(c - c')/i(c) 



(10) 



where, in the last step, we introduced the static pair correlation function g by /2(r — r', c, c') = fi(c) fi( c')g (r — r'). 
As the static structure factor S(r) = (Sp(r + r)Sp(r Q )} is related to the pair correlation function by |36l. |37| 

5(r) - P l(g{v) - 1) + w 5(r) , 

the desired correlation function for the fluctuations of the one-particle distribution function is finally obtained as 

(5/(r,c)5/(r',c')) = [/»/>')[%- r')/p - ^(r - r')}/p + ^/(c)*(c - c')5(r - r')] . (11) 

The first term on the r.h.s. of (fTTj) describes spatial correlations due to the non-ideal character of the fluid. For 
the ideal gas, we have g(r) = 1, S(r) = ppo5(r), thus only the last term remains, and relation (|TT|) becomes the 
well-known expression used in the Boltzmann-Langevin theory for a dilute gas [2^, [3(J HH , 

(<5/(r, c)Sf(v', c')> = pf(c)S(c - c')5(r - r') . (12) 
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The presence of the two delta-functions reflects the fact that in the ideal gas, there are no correlations between two 
particles at different positions or with different momenta. The factor p arises due to the definition of the distribution 
function in terms of the mass density. 

It can be checked that expression (Tnj) is self-consistent regarding the definition of the structure factor: As the 
density is the zeroth moment of the distribution function, p(r) = J f(r,c)d d c, we immediately obtain (Sp(r)Sp(r')) = 
J d d cd d c' (6f(r, c)Sf(r', c')) = S(r — r'). Similarly, by the definition of the macroscopic fluid momentum as the 
first moment of the distribution function, jo,(r) = J d d cc a f(r 7 c), we obtain the well-known expression Q for the 
equal-time momentum correlation function, 

(6j a (r)6 Jp (r')) = J d d cd d c'c a c' (Sf(v,c)Sf(v',r')) = p, J d d cc 2 J(c)5(r - v')5 afj = Po k B TS(r - r')5 af3 , (13) 

Here, we used the fact that / [eq. 10 ] describes a quiescent state, i.e. J d d ccf(c) = 0, and the Gaussian integral 
J d d cf(c)c a cp = SappokfiT / p. Expression (|13[) shows that, both in the ideal and non-ideal fluid, there are no 
correlations in the momenta of different fluid elements. Indeed, as is well-known from statistical mechanics, the 
Hamiltonian of an equilibrium fluid is diagonal in the momenta. 



III. DISCRETE VELOCITY BOLTZMANN EQUATION 



In the derivation of the lattice Boltzmann equation, the kinetic space is discretized in the velocities to yield a 
discrete velocity Boltzmann equation for the distribution function f(r,c,t) (39j . Choosing an appropriate set of 
discrete velocities Cj, the discrete velocity Boltzmann equation for the distribution function /j(r, t) = f(r, c,, t) under 
the action of an arbitrary body force F reads 

dtfi + c, • V/i + [F • V c /]i = -An {fj - f- q ) ■ (14) 

Here, fi(r,t) represents the mean number of particles at r and time t, moving along the lattice direction defined by 
the discrete velocity a, = 1, . . . , N). Furthermore, 

,cq / . - Ci pv a vpQ ia0 \ 

fi =Wi { P+ ^r + 2 of ) (15) 

is a local equilibrium distribution, the discrete analogue of a local Maxwellian distribution in continuum kinetic 
theory truncated to second order in the mean flow velocity v, the Wi are a set of weights which satisfy J^- u>j = 1, 
and Qi a p = Ci a Cip — cr^S a p. Greek indices denote Cartesian directions and summation over repeated free indices is 
implied. The constant a s is the speed of sound of an ideal LB gas, and its numerical value is fixed by the structure 
of the discrete velocity set. For the conventional D2Q9 lattice employed below, it has a numerical value of a 2 s = 1/3 
in lattice units. Note that this quantity has to be distinguished from the actual speed of sound c s of the fluid under 
consideration. If the fluid is non-ideal, c s will in general be different from a s . The quantity [F • V c /]i represents the 
discretized equivalent of a forcing term, which we take as (ioj 

iiS - [P .VJ] i =p» i (^ + (™^i). (16) 

For present purposes, where a body force is supposed to mediate the non-ideal fluid interactions in the fluctuating 
DBE, we will assume that F depends only on the density and its derivatives, F = F(p, Vp). 

The low order velocity moments of the distribution function are related to the densities of mass, momentum and 
the deviatoric stress, 

{p,pV a , S a p} = } J fi{h Cj a , Qj a p} (17) 
i 

where S a p + pa 2 S a p = Tl a p is the Eulcrian momentum flux. The higher moments of the distribution function are 
related to the densities of rapidly relaxing kinetic degrees of freedom, variously called ghost or kinetic variables. 
In eq. (fl"4|) . Ay is a scattering matrix whose eigenvalues control the relaxation of the kinetic modes to their local 
equilibrium values [H, [l]| . The null eigenvalues correspond to the eigenvectors associated with the conserved mass 
and momentum densities, while the leading non-zero eigenvalue associated with Qi a p controls the viscosity of the LB 
fluid. 
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For a general athermal DBE model with n velocities in d space dimensions, the n x n collision matrix has d+ 1 
null eigenvectors corresponding to the conserved density and d components of the conserved momentum, d(d + l)/2 
eigenvectors corresponding to the stress modes, and n — (d + 1) — d(d + l)/2 eigenvectors corresponding to the 
ghost modes [HI, |4l[. The choice of the null and stress eigenvectors {1, Cj Q , Qiap} follows directly from the physical 
definition of the densities associated with them. Without specifying the exact analytical expression for the remaining 
eigenvectors, let us assume a linearly independent set of the eigenvectors of the scattering matrix be given by {T a i}, 
where a = 1 . . . n labels the eigenvector, and i = 1 . . . n labels the component of the eigenvector along the ith velocity 
direction. This allows us to define densities associated with the eigenvectors T a as moments of the populations by 

m a (r,t)=T ai fi(r,t). (18) 

For T ai € {1, Ci a , Qia/s}, the densities are the mass, momentum and stress. The ghost eigenvectors are higher 
polynomials of the discrete velocities. The discreteness of the kinetic space implies that, unlike in the continuum, only 
a finite number of polynomials can be linearly independent, being equal to the number of discrete velocities. For a 
model with n discrete velocities, the choice of the n linearly independent polynomials is not unique, but defined only 
up to a similarity transformation. Independent of the precise choice, the distribution function itself can be expanded 
in a linearly independent set of eigenvectors which are polynomials of the discrete velocities 

fi(v,t) = Wi^m a (r,t) . (19) 

Consistency between the above two equations implies that the set of polynomials T a arc both orthogonal and complete, 

J2wiT ai T bi = N a S ab , (20) 



£ ^ = . ( 21 ) 



where N a is the length of the ath eigenvector. Crucially, with the definitions above, the eigenvectors T a i form an 
orthogonal set under a weighted inner product (T a ,T D ) = Yli w iTaiTbi [H|- Compared to using an unweighted inner 
product [il|, , the advantage of the present choice is that the equilibrium distribution (fl~5|) has no projection onto 
the ghost sector, i.e. T a if° q = if T a is a ghost eigenvector. This emphasizes the physical interpretation of the 
model, but is otherwise not necessary for a faithful implementation of hydrodynamics [431 ] . It is often convenient to 
imagine the fi to be the elements of vector in the space spanned by the velocity vectors Cj . Then, T a i defines via (fT8|) 
an orthogonal transformation matrix between the distribution function-space and the moment-space. Relation (|19[) 
serves to define the transformation matrix in the inverse direction, (T _1 ) !a = WiT a i/N a . One possible choice for the 
T a for a D2Q9 lattice is presented in appendix [Bj 

Based on the notion that T a represents polynomials corresponding to the hydrodynamic, transport, and ghost 
eigenvectors, correspondingly, eq. (|19[) allows the distribution function to be separated into contributions from the 
hydrodynamic, transport, and ghost moments 

fi =f? + fi+fi- (22) 



This structure can be conveniently organized by introducing projection operators [4J] which project the distribution 
function onto the hydrodynamic, transport, and ghost subspaces, 

~2^f 



P?fj = fi = W . (24) 



N a 



The explicit form of the projection operators immediately follow from the completeness relation (|2lj) as 

pH _ ( -, , C '' C j 1 pT _ w iQiupQjaP pG _ \ " w iT a iT a j 

aeG 



2(J 4 > U ^ N a 



The discrete Maxwellian and the forcing term is a nonlinear function of the distribution function, and thus cq. (|14[) 
is only apparently linear, the nonlincarity being concealed in /° q and $i. A useful linearization of the LB equation 
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consists of neglecting the quadratic term in the discrete Maxwcllian to yield a local equilibrium h eq which is linear in 
the mean velocity, 

C = + =if/i- (26) 

Similarly, the forcing term $j, eq. (|16p . can be linearized as 

<fH = Pfi$j(SF) = pwi^i = . (27) 

where the last equality serves to define the linear forcing operator K^j. This definition is always possible, as the 
linearized interaction force SF is assumed to depend linearly on the density and its derivatives. While the definitions 
of SF and are formal at the level of eq. (f2~T)) , their precise forms will become clear in the following section, where we 
discuss a particular DBE model for a non-ideal fluid. In the linearized approximation for the equilibrium distribution, 
we have f^ q = h eq = (P H /),;. Defining a right-projected relaxation matrix as Kf- = Aifc(l — P H )kj, the relaxation 
term of the DBE (fH| can be written as 

M/i-/r]=A« [h-{P H f) j \=Kf j f j , 

which is obviously linear in Thus, the linearized DBE follows as 

BtSfi + a ■ VSfi = -AfiSfj + KySfj (28) 

where the obvious notation Sf% is introduced to indicate a fluctuation. 

It is clear that by construction has eigenvectors of mass and momentum with zero eigenvalues. The form of the 
matrix, by itself, places no constraint on the eigenvalues of the transport and ghost sectors. The simplest possible 
choice consists of the BGK approximation, where the relaxation matrix is diagonal, Ay = Sij/r, which implies that all 
non-conserved modes relax at the same rate 1/r. In general, the maximum number of different eigenvalues depends 
on the dimensionality and the chosen basis set T a The fluctuating DBE that we will present below places no 
further constraints on the eigenvalues, and we will denote the eigenvalues of A^ simply by A a , i.e. 

A%T aj = X a T ai , (29) 

keeping in mind that by construction A a = for the conserved modes (i.e. a = 1, . . . , d + 1). 

The linearized DBE (|28[) is most conveniently analysed in Fourier space. Introducing the Fourier transformation 
by 6fi(r) = (I/27r) d / 2 / Sf^k) exp(-ik ■ r) dk yields 

d t Sfi(k) - ik ■ ci*/<(k) = [-A£ + Ay(k)] 6f 3 (k) . (30) 

Thus, the DBE dynamics is reduced to a set of coupled ordinary differential equations for the Fourier modes Sfiik, t). 
The linearized DBE can equally well be written in terms of the Fourier transforms of the moments Sm a (k,t), using 
definition (fTSl) to obtain 



d t 5m a (k) + A ab {k)8m b {k) = -\ a 5 ab + K ab {k) Sm b {k) , (31) 

where A ab (k) is the Fourier-transformed advection operator (see appendix [C]) 

A ab {k) = -ik • Ta.c^- 1 = -ik • WiT ai T biCi /N a , (32) 

and K ab = TaiKijT^ 1 denotes the forcing operator in moment space. Clearly, since the populations and densities are 
related by an invertiblc linear transformation, eq. (|3I[) represents an exact reformulation of the DBE dynamics in the 
basis of moments. The dynamical equation in the /, basis diagonalizcs the advection operator — ik ■ Cj&y, while the 
dynamical equation in the m a basis diagonalizes the collision operator A^ . The eigenvectors of the dynamics are a 
combination of the /j and the m a . A forcing term will contribute off-diagonal terms in either representation, as the 
m a are constructed as eigenvectors of the relaxational part of the collision operator, A^. 
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IV. FLUCTUATING DISCRETE- VELOCITY BOLTZMANN EQUATION 

We shall now turn to the central result of the present work: the fluctuating DBE and the corresponding fluctuation- 
dissipation relation. The fluctuating DBE is obtained by promoting the DBE (|2"5)l into a Langevin equation, where 
the fi are interpreted as instantaneous, fluctuating densities in phase space (corresponding to the // in the notation 
of section HTCj) . 

dt.fi + Ci ■ Vfr = -Aytfj - jf ) + <!>,+ i, . 

The noise terms £j are Gaussian random variables that give rise to fluctuations in the populations in each phase space 
cell. Restricting to linear fluctuations, the FDBE becomes 

dtSfi + ^ ■ VSfi = -A* Sfj + Kij 5fj + & , (33) 



where, by construction, the linear relaxation matrix Kf^ has zero eigenvalues for the conserved modes [see eq. (|29|) ]. 
After Fourier-transforming the linear FDBE (|33j) and expressing it in terms of moments, we obtain 

d t 5m a {k) = \-A ab {k) - \ a 6 ab + K ab (k)] 5m b (k) + f a (k) 

L _ J (34) 

= -L ab (k)Sm b (k) + £ a (k), 
with £ a = Tai^i. For later convenience, we introduced the DBE time-evolution operator L by 

L ab {k) = A ab (k) + X a 5 ab - K ab (k) . 
Above linear FDBE can be brought into the form of Langevin eq. © by writing it as 

d t Sm a (k) = - j dk'£ ah (k, k')m 6 (k') + £ a (k) , 

with C ab (k, k') = L a f,(k)«5(k — k'), and interpreting the wavenumber arguments formally as indices. Hence, the 
operator defined by X a 8 ab S(k — k') is Hcrmitian in the wavenumber and mode indices, while the operator defined by 
A ab (k)S(k — k') is anti-Hcrmitian in the wavenumber indices [cf. (|32|) ]. but of mixed character in the mode indices. 
As will be shown below, the same holds for the operator K a b{k)S(k— k'). Denoting the equal-time correlation matrix 
of the modes by Q ab (k,k') = (rn a (k)m b (k')) , the noise correlations follow straightforwardly from the FDT, eq. ((H), as 

(Uk,t)i b (k',t')) = J dk" [£(k 1 k")G(k",k') + g(k,k")tf(k",k')] ab 6(t-t'). (35) 

The above relation can be simplified by noting that, due to translational invariance, the mode correlations are generally 
diagonal in Fourier space, i.e. Q(k, k') = G(k)S(k + k') with G ab (k) = (m a (k)m b (— k)) = (m (k)m^(k)). This allows 
to obtain the FDT in its final form 

S ab (k) ee (i(k)| 6 (-k)) = [L(k)G(k) + G(k)Lt(k)] afe , (36) 

where we defined the equal-time correlation matrix of the noise S ab (k). Within our convention, Gn is the equal- 
time density correlation function (static structure factor) and G aa for a = 2,...,l + dis the equal-time momentum 
correlation function. 

In contrast to the time-evolution operator L, which is provided by the DBE, the correlation matrix G must in general 
instead be derived from a statistical mechanical framework. The specification of G constitutes a central ingredient of 
any FDBE model. While the matrix elements of the conserved subspace of G can immediately be determined from the 
basic statistical mechanical theory of fluids Q, the specification of the remaining elements is less straightforward. In 
particular, it must be emphasized that the correlations of the non-conserved modes are not arbitrary, as the full matrix 
G must combine in the FDT (f31))) with the reversible and irreversible parts of the DBE dynamics (as represented 
by L) to obtain a positive-semidefinite noise covariance matrix for all wavenumbers. The intimate relation between 
equilibrium correlations and dynamical equations can be seen from the following argument based on the FDT of the 
fluctuating Navier-Stokes equations of a non- ideal fluid [H EH, EH : due to the coupling between density and momentum 
fluctuations, once the expression for the random stress tensor correlations has been found, the structure factor can 
be determined from the dynamical equations without any further reference to statistical mechanics. This illustrates 
the general fact that driving forces in the dynamical equations arc provided by thermodynamic forces and therefore, 
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already contain information on the equilibrium fluctuations [23|, |33[. Transferred to the DBE, this implies that the 
expression for the structure factor is related to the interaction force and, more generally, the equilibrium correlations of 
the conserved modes are indirectly fixed by the correlations of the non-conserved modes due to the dynamic coupling 
between both. As shown in [l6| in the context of the LBE, in principle, it suffices to derive a fluctuation model for the 
non-conserved modes and obtain an expression for the noise based only on this information using a detailed balance 
argument. In the context of the Langevin approach to the DBE, however, it will be necessary to specify the full 
equal-time correlation matrix of all modes, paying special attention to the relation between equilibrium correlations 
and interaction forces. 

Finally, it should be noted that the form of the interaction force is in principle restricted by the numerics of the 
streaming process. This can be seen, for instance, by considering a situation of coexisting phases in a quiescent 
equilibrium fluid: in order to maintain the Maxwellian equilibrium distribution (with spatially varying density) at 
each point in space, the non-zero contribution of the forcing term in the collision step of the Boltzmann equation must 
be compensated solely by the streaming term, as the collision term identically vanishes. This requirement becomes 
delicate to fulfill in the LBE where space is discretized and demands a careful implementation of the forcing term. 

In this work, we adopt eq. f|l If) . which is directly based on continuum kinetic theory, as the defining relation for 
the equal-time correlations of the distribution function in both the ideal and non-ideal fluid. Written in terms of the 
discrete velocities, eq. (|11[) becomes in Fourier-space 

(*/i(k)*/j (k')> = [fjj [S(k)/ Po - M ]/p + Ufa] S(k + k') , (37) 

where fi(c) = f^ q (c, po,v = 0) is the discrete Maxwellian. Due to translational invariance, all correlation functions 
are diagonal in Fourier-space, thus allowing us to write them in terms of a single wavevector. The moment-space 
correlation matrix is obtained by projecting (|37[) onto the basis vectors, 

G ab (k) = (5m a (k)Sm b (-k)) = T ai T bj (Sfi(k)Sfj (-k)) 
= (m a m b [S(k)/p - (j]/p + pT ai T bl ji) , 

where Tn a — T a ifi. Whenever the equilibrium distribution is of ideal gas form, ec{. (|15p , we have fn a — po^ai* 

In the 

above relations we have used the fact that for the Fourier-transform of a real- valued quantity, such as the distribution 
function or its moment, complex conjugation is equivalent to changing the sign of the wavevector. 

The structure factor S(k) in general depends on the underlying (non-ideal fluid) interactions and must itself be 
determined from a statistical mechanical theory, i.e. it is not provided by eq. (f3"T|). Note however, that for the force- 
based DBE to be discussed below (section IIV B 1[) , the structure factor will drop out of the final expression for the 
noise covariance S due to the reversible nature of the interaction force. In contrast, this will not be the case for the 
modificd-cquilibrium model (see section HV B 2[) . 

Crucially, in contrast to the continuum Boltzmann case, the parameter p representing the mass of the fluid particles 
is not known anymore a priori in the algorithm, but now must instead be determined self-consistcntly from the 
expression for the momentum correlation function obtained from eq. (|37p . The discrete momentum correlation function 
of any fluid in equilibrium at a temperature T and density po is assumed to be given by the same relation (| X3[) as in 
the continuum, 

(6j a (k)Sjf> (k')> = p k B TS aP S(k + k') , (39) 

where the momentum density is defined as j = /jCj [see eq. (|17[) ]. From this relation, we will find that p is related 
to the fluctuation temperature of the fluid, in agreement with a previous analysis of the fluctuating ideal gas LBE 
[l6|. This is explicitly shown below (sees. ITVAl and TlVBp . where we discuss the ideal and non-ideal fluid models. 
Thus, once the structure factor and the equilibrium distribution are specified, relation (|38p allows to compute any 
correlation function of two modes. 

In the FDT, eq. (|3l)|) . it is assumed that the noise is not correlated in time, i.e. the £ a represent Gaussian white noise 
in time. This is physically well justified, as the noise represents the "fast" degrees of freedom of the many-particle 
system which typically have correlation times many orders of magnitude smaller than the interesting macroscopic 
quantities @,|46|. Similarly, one might reason that the noise should also not be correlated in space, as the fast variables 
relax on scales much smaller than inherent to the description of the Boltzmann equation. Spatially uncorrelated noise 
will indeed be an exact result of the FDT for the force-based models discussed below. The model of Swift ct al. 
[l9l ]. however, deviates from this basic intuition in that it requires spatially correlated noise, as shown in [l8| . The 
origin of this result can be traced back to the use of a modified equilibrium distribution, which is also responsible for 
the non-local (wavelength-dependent) bulk viscosity arising in the associated Navier-Stokes equations of that model. 
Besides the absence of temporal correlations, no further constraints will be imposed on the noise. In particular, the 
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laws of mass and momentum conservation are not enforced a priori, but will instead be an outcome of the FDT by 
construction of the relaxation matrix A R , as shown below. 

Finally, we shall make some comments on the applicability of the noise derived from a linear theory to situations 
where non-linear terms are present. As can be shown on general grounds via a Fokker-Planck treatment [36l I45ll47|. 
in cases where a well-defined equilibrium distribution (e.g. in terms of a free-energy functional) exists, both the noise 
obtained from the linear and fully non-linear hydrodynamical treatment are identical. Thus, the noise expressions 
derived for the models considered in this work are expected to remain valid in the presence of non-linear couplings, 
provided they are either fully reversible or stem from an underlying free-energy functional. The advective non-linearity 
jajp/p of hydrodynamics belongs to the first class, while cubic (or higher order) density terms in the interaction force 
(equation of state) represent so-called dissipative couplings originating from a free-energy functional Q . Although a 
strict treatment should deal with the appropriate Fokker-Planck equation for the FDBE, we nevertheless expect that 
the noise derived in the present work can be used to model fluctuating non-linear fluids, as, for example, required for 
phase coexistence or critical phenomena. 

Under general non- equilibrium conditions, as, for example, a fluid under flow, the situation is not so straightforward. 
Here, it has been shown that the equations of fluctuating hydrodynamics remain valid in simple non-equilibrium steady 
states (e.g., in a temperature gradient or uniform shear), provided that the noise covariances are computed with the 
local values for the thermodynamic variables and transport coefficients 0, [46|, I48l - l53j . This is intuitively clear from 
the fact that the fast variables are uncorrelated on macroscopic length and time scales and can thus be expected to 
be in a local equilibrium. Thus, it seems reasonable to assume that FDBE noise with a covariance of the same form 
as in the equilibrium case can also be employed in simple non-equilibrium states. However, we are not aware of any 
studies investigating these issues in detail for the FDBE or FLBE at the present time. We remark that violation of 
Galilean invariancc appears to limit the application of the FLBE to systems with rather small mean flow velocities 
[54| . In this work, we will be only concerned with fluctuations around equilibrium states. 



A. Ideal Gas 



As a first application of the theory presented above, we derive the noise covariance for the ideal gas model. Our 
result, obtained for the FDBE, will be shown in section [V] to be identical to the result derived by [l5[ using the direct 
FLBE approach. The FDBE of the ideal gas is given by 

d t 5m a (k) + A ab (k)8m b (k) = -X a Sm a (k) + £ a (k) , (40) 

and, thus, the generalized time-evolution operator of eq. Q34p consists in this case only of an advective and a relaxational 
part, 

L ab (k) = A ab (k) + X a 6 ab , (41) 

where A a= i a+i = by construction of the relaxation operator. To compute the equilibrium correlation matrix of 

the modes G, we assume the structure factor to be given by the conventional (wavelength-independent) ideal gas 
expression S'id(k) = p$p, with a yet undetermined mass parameter p. Hence, we obtain from eq. (|57|) the DBE analog 
of relation (|12p as 

(SM^Sf^k')) = pf t 6 tJ 6(k + k') , (42) 

where fi = poWi is the global Maxwellian. Clearly, the independence of (|42|) of the wavenumber expresses the absence 
of all correlations in an ideal gas. With the help of the orthogonality relation (f2"U)) , we compute the equilibrium 
correlation matrix as 

G ab (k) = (<5m a (k)<5m & (-k)) = T ai T b j(SfiQ<.)5fj(-'k)) = ppaw l T ai T bi = pp a N a 5 ab , (43) 

where N a is the length of the ath basis vector. The parameter p can now be determined by requiring consistency 
with the basic statistical mechanical relation (|39[) for the momentum correlation, fixing 

M=^, (44) 

since N2,...,d+i = a 2 s (c.f. Tablc|l|. As a s is also the speed of sound of the ideal gas, we see that (|44|) is consistent with 
the ideal gas equation of state if p is interpreted as the mass of a fictitious DBE-particle [l6|. From (|43|) together 
with the above value of p the complete expression for the structure factor of the ideal gas is obtained as 



(45) 
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With relations (|4T1) and (|43j) . the noise covariancc matrix finally follows from eq. (|36j) as 



E ab (k) = 2X a G ab (k) = 2\ a 



p k B T 



N a 5 ab 



(46) 



Transforming back into real-space and invoking the definition of S a b(k) as (£ a (k)£f,(k')) = H a b(k)5(k + k'), we obtain 

(6,(r)&(r , )> = 2\ a ^^NJ ab 5{v - r') , 

which explicitly indicates the absence of spatial correlations. Crucially, as A Q = by construction for the conserved 
modes [see eq. f|29|>]. noise with the above covariance automatically obeys mass and momentum conservation, i.e. 
£ a = for o=l,...,d+l. Notably, the equilibrium correlation matrix of the ideal gas and the advection operator 
fulfil AG = -GA*, which is the technical reason why the above noise covariance is independent of any advective 
contribution, in line with physical expectation. This was also found in (l5j using the FLBE approach to the ideal gas. 



B. Non-Ideal Gas 



We now proceed to analyse thermal noise in the DBE of two non-ideal gas models. First, as a central new result, 
a general force-based FDBE is presented that applies, among others, to the model introduced by He, Shan and 
Doolcn 24 1 . Next, the modified-cquilibrium model of Swift et al. [19|,[55| is analysed within the DBE approach, thus 
complementing previous work based on the LBE [l8j . The equations of fluctuating hydrodynamics of non-ideal fluids 
resulting from the two FDBEs below can be found in |l8| . 



Force-based model 



The general FDBE for a force-based non-ideal fluid model is given as 

<9 t <5m a (k) + A ab (k)Sm b (k) = -\ a 5 ab + K ab (k) 5m b (k) + £ a (k) 



(47) 



where K ab is the linearized forcing-term defined through K ab m b = T a ifa with fa given by eq. (|27p. fa = poWi6F lnt -c/ag. 
The equilibrium distribution is assumed to have the standard Maxwellian form, eq. (|15[) . In the linear approximation, 
the forcing-term only affects the momentum density, since 



T ai fa = {0,6F a ,0,...}. 



(48) 



To proceed, an expression for the interaction force in terms of the density has to be specified. We shall assume the 
linear interaction force to be given by 

5F^ = ik[c? s (k)-^]6p, (49) 
where c s (k) is a generalized speed of sound that is related to the structure factor of the non-ideal fluid by 

Pok B T 



S(k) 



(50) 



It is shown in appendix [Al that assumptions (|49|) and ((50)) hold, in particular, for the model of He, Shan and Doolen 
[Hj|. Crucially, for the purpose of deriving the FDT, the specific expressions of c s (k) or S(k) are not important. 
Hence, the present model is not restricted to a particular thermodynamic framework, such as a square-gradient free 
energy functional, used to describe density fluctuations. 

Combining (|48|) with (|49|) . we see that the forcing operator in the FDBE (|47|) of a general D2Qn model can be 
written as 



K ab (k) 



iky 




"c2(k)-*f 
^(k)-a s 2 " 







\ 
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where the upper- left submatrix spans over the conserved modes. The generalization to the three-dimensional case 
is obvious. Using the fact that the equilibrium distribution is of ideal gas form, fh a = Pq8 ci i, the mode correlation 
matrix G follows from eq. (f3"8j) as 

G(k) = diag [S(k), ■ ■ ■ , WoN n ] , (51) 

with the structure factor being given by eq. (|50[) . From eq. (|51[) we immediately see, that the mass parameter p, has 
to be chosen as 



/' 



in order to obtain the correct momentum correlations, G22 = ■ ■ ■ = Gd+i.d+i = Po^bT- Thus, p is identical in the 
ideal and the force-based non-ideal fluid model. This is not surprising, as both models employ the same equilibrium 
distribution. With the time-evolution operator being given by 

L ab (k) = A ab (k) + X a S ab - K ab (k) , (52) 

the noise covariance matrix follows from the FDT (|36| after a bit of algebra as 

~ ab (k) = 2A a G Qb (k) = 2\ a ^^-N a 5 ab , (53) 

where by construction A a = for the conserved modes (a = 1, . . . d + 1). 

The above noise covariance is identical to the one obtained in the ideal gas model, eq. (|46| , implying that the non- 
ideal interactions have no effect on the dissipation. This is physically expected as the forcing should represent a fully 
reversible contribution to the DBE dynamics. Technically, this result depends crucially on the cancellation between 
terms originating from the interaction force by corresponding terms originating from the linear advection operator in 
the FDT. An fundamental prerequisite for these cancellations to occur is the assumption that the interaction force 
is of the general form given by eq. (|49|) . with the speed of sound being related to the structure factor by (|50| . If 
these prerequisites are met, both contributions of forcing and advection to the FDT, eq. (f3"r?)) . disappear from the 
final result, eq. (|53[) . Conversely, the requirement to obtain a well-defined FDT obviously also allows one to impose 
constraints on the possible forms of the interaction force. 



2. Modified- equilibrium model 

The Langevin extension of the modified-equilibrium model of Swift et al. [ill [55| has been analysed in [l8| based on 
the LBE dynamics. There, it was shown that relation (|37p for the equal-time correlations of the distribution function 
had to be modified to obtain a well-defined noise covariance matrix. The latter turned out to be fc-dependent and 
contained residual influences of the thermodynamic interaction model. It will be informative to analyze the modified- 
equilibrium model starting from the DBE, and thereby compare it to the force-based model of the previous section. 
Where necessary, we will specify expressions in the D2Q9 basis given in appendix [B] for simplicity. 

To derive the FDBE of the modified-equilibrium model, we note that in this case the equilibrium distribution / eq 
is not just given by the projection of the full distribution onto the hydrodynamic subspace, P^ fi [eq. (|26|)]. as the 
second moment of / oq is defined to reproduce an equilibrium non-ideal pressure tensor. The derivation leading to 
eq. (EHD can, however, easily adapted by noting that A(/ - f cq ) = A(/ - P H f) + A{P H f - f cq ) = A R f - A R f cq , since 
P H f eq = P H f due to mass and momentum conservation. Thus, in moment-space, we arrive at the FDBE 

d t Sm a (k) + A ab (k)5m b (k) = -X a 6 ab [Sm b (k) - 5m cq (k)} + | (k) , (54) 

where again, A a = for the conserved moments and Sm a q is the linearized equilibrium distribution in moment space. 
In this model, all non-ideal fluid interactions are represented by a pressure contribution to the e quili brium distribution. 
Specifically, for the D2Q9 basis set defined in appendix [Bl the equilibrium moments follow as [18[ 

Sm cq (k) = {5 Pl Sj x , Sj y ,d(k)Sp, 0, 0, 0, 0, -d(k)5p} , (55) 

where d(k) = 6 [c^(k) — af\ and c s (k) is a generalized speed of sound which is related to the structure factor S(k) 
by the eq. (|50|) above. As seen in expression (|55|). the non-ideal equilibrium contributions appear in a bulk pressure 
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and in a ghost mode. Similar results are obtained also on three-dimensional lattices. It was shown in [ljj], that the 
modificd-cquilibrium model requires the following equilibrium correlation matrix: 

(5f i (k)Sf j (-k)) = ^-f i (k)6 ij . (56) 
Po 

This is different from the expression cq. (|37[) derived from continuum kinetic theory. The reason for these discrepancies 
can be traced back to the fact that, in contrast to conventional kinetic theory, the modificd-cquilibrium model employs 
a non-Maxwellian form of the equilibrium distribution. With cq. (|55[) . the time-evolution operator of the DBE (|54|) 
can be identified as 

L a b(X) = A ab (k) + X a S ab - \4,d(k)5 a4 ,6bi + Xgd(k)S a9 S b i . (57) 
The noise covariance matrix then follows from the FDT ([36]) after a bit of algebra as 
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U[c 2 s (k)-a*}(A e + \ e ) . 


N 7 X q . 
■ N 8 X q 

iV 9 [|_| c 2 (k)]Ae J 



(58) 



where the dots indicate zeros for short. 

While the noise obtained for the modificd-cquilibrium model exactly respects mass and momentum conservation, 
it is wavcnumbcr-dcpendcnt and, hence, non-local, in contrast to the noise (|53[) of the force-based model. These 
non-localities arise due to the modified equilibrium distribution employed in the model, c.f. cq (155|) . which is also 
responsible for the wavenumber-dependent bulk viscosity arising in the hydrodynamic equatio ns UM . Expression (f58|) 
is identical to the noise obtained from the LBE analysis of the modified-equilibrium model jl8| . except for lattice- 
induced contributions to the relaxation parameters (see below). Similarly as for the ideal gas, we find that also 
for the modified-equilibrium model, the equilibrium correlation matrix and the advection operator fulfil the relation 
AG = -GAt Hence, the contribution of the advection operator to the FDT identically vanishes, in contrast to the 
force-based model, where the advective contribution is canceled by a contribution from the forcing term. 



V. FLUCTUATING LATTICE BOLTZMANN EQUATION 



In order to transfer results of the previous section to the LBE, we apply a second-order accurate characteristics 
based integration to the FDBE, following J5||. The subsequent steps can be performed on the general (non-linear) 
FDBE 

dtfi + c, ■ Wfi = -Aij (fj - /f) + (59) 
which we write in compact notation as 

d t f i + c i -Vf i = R i (r,t), 

where Ri — — Ay-(/j — /| q ) + $i + & is introduced for short. Integrating over a time step At, we obtain 

rAi 



fi (r + Ci At, t + At) - f. t (r, t) = dsR, (r + as, t + s) 

Jo 



= ^Ri(r + ctAt, t + At)- -^Ri{r, t) + AtR t (r, t) . 

Evaluation of the integral using the trapezium rule leads to a set of implicit finite-difference equations for the fi, 
which can be made explicit by introducing a set of auxiliary distribution functions 

ft B (r,t)^f t (T,t)-^-R t (v,t). (60) 
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The time step is set henceforth to At = 1 in lattice units. Expressing Ri in terms of the / 2 LB , one obtains 

-l 



Ri = 



which leads to the fluctuating LBE (FLBE) 



-K jk (f k 



LB 



•J 



fl B {v + Cll t + l) = fl*{v,t)+(l + -k) [-KM 



f LB 
Ik 



(61) 



(62) 



Note that the noise variables appear here just as another source term in addition to the forcing-term In order 
to write the LBE in moment space using an arbitrary basis set T a , we compute (I + ^A) _1 = fl + ^T~ 1 AT 

T, where in the last equation, we defined A = diag[A Q ], A LB = diag[A^ B ] and 

LB 



T 



T = T~ 



introduced a set of collision parameters A„ by 

A^ B = 



1 

Xb 



1 



2X r , 



T a + 



2 + A Q 



(63) 



These redefined collision parameters allow us to write the LBE in a form conventionally found in the literature. With 



these definitions, we obtain (1+ iA) _1 A = T _1 diag[2A a / (2 + A a )]T = T~ 1 K LB T, leading finally to the moment-space 



version of the FLBE (|62|). 

fl B (r + Ci ,t+l)=T^[m^- 



LB 



LB 



1 - -A 



LB 



(64) 



We recognize in (|64|) the appearance of the well-known factor (l — ^Aq B ) in front of the forcing-term |56l - [53 |. which 
is necessary to ensure a second-order accurate influence of the body force. Since our derivation of the FLBE does 
not differentiate whether source terms originate from the random noise or a body force, the same factor naturally 
multiplies also £ Q . Using eqs. ([60]) and (fSTj) . the original DBE moments m a = T a [fi can be expressed in terms of the 
redefined LBE moments m^ B = T ai fl B by 



LB 



1 



-Ar Kr - < q ) + ( 1 - 2 x a ) « + ^ 



In particular, the hydrodynamically relevant moments density, momentum and stress are obtained as 

fi > 



LB, 



Sap = f 



LB, 



(iap 



2 

1 

2 



-K (ft "Qi*P - P^aUp) +[l--K 



LB 



(u a Fp +UpF a + ZiQiap) 



(65) 

(66) 

(67) 
(68) 



where in the last equation, the expression for the second moment of the forcing term, eq. (|16|) . has been used. Note 
that due to eq. ([53]) . the noise gives a non-zero contribution to the instantaneous stress S a p- The last two equations 
contain the well-known lattice corrections to the momentum and stress in force- based models p56T - [58j . 

The term (1 + \\\ B )i a = ia B in eq. (JUJ) defines the noise |^ B appropriate to the FLBE in terms of the FDBE 
noise £ . Both the ideal and the force-based non-ideal fluid FDBE considered above have identical noise covariances, 
given by eqs. ([46]) and (|53]l. Hence, in both cases the real-space covariance of the LBE noise follows as 



1 



(C LB (r)^ B (r')> = 1 - TiKl 1 - ^ B (Ur)W)) = 



1 



ol AV 



(2 - \™)\™ N a 5 ab 5 T 



(69) 



The factor AV arises from the lattice equivalent of the delta function and is taken as the volume of the elementary 
lattice cell. It reflects the fact that fluctuations become more pronounced with decreasing length scale. Expression 
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(|69|) is identical to the results previously derived by jTH], [T^| for the ideal gas LBE 2 . Comparing eq. ([69]) to (|46|) . we 
identify the term — (A^ B ) 2 as a lattice correction to the FDT of the FLBE, analogous to the well-known "streaming 
contribution" to the viscosity. 

In contrast to the modified-equilibrium model, the derivation of the FDT for force-based non-ideal fluid models 
turned out to be technically cumbersome using the direct LBE approach introduced in (l8j . There, difficulties arise 
from the presence of the non-linear LBE advection operator, A LB (k) a b = T a j exp(— ik • Cj)T^ , and the spatial 
discretization scheme of the interaction force. In contrast to the DBE, it appears to be technically demanding to 
construct a LBE forcing term that is fully reversible at the lattice level, while at the same time fulfil basic physical 
requirements (such as thermodynamic consistency (60l. l6lj or the absence of spurious currents [(32|). In this reg ard, 
the DBE-based approach, as presented in this work, cannot fully replace a lattice treatment of the FDT |15l - [l8j . 
Nevertheless, it provides a useful starting point for the construction of fluctuating LBEs. 



VI. SIMULATIONS 



We now investigate whether the fluctuating LBE of the force-based model derived in the preceding section can 
correctly reproduce some basic statistical mechanical results for non-ideal fluids. Results for the ideal gas and the 
modified-equilibrium model have already been reported in [l5j and [l8| . where close agreement between simulation 
and theoretical expectations has been found. Here, we first check whether thermal noise defined by eq. (p9"|) leads 
to the correct equilibration of all degrees of freedom in a LB simulation of a homogenous non-ideal fluid. Next, 
capillary fluctuations in a liquid-vapor interface are investigated as an important example for thermal fluctuations in 
an inhomogeneous system. For the generic force-based model discussed in this work, spatially uncorrelated, ideal-gas- 
like noise with variance given by cq. (|69|) is an exact consequence of the FDT for all wavenumbers. This form of noise 
is easily implemented in a simulation, as all noise modes £ a can be drawn independently on each lattice site from 
a Gaussian distribution. For this purpose, the fast random number generator described in [63| can be successfully 
employed. 

Simulations of the non-ideal fluid are performed using D2Q9 implementation of the model of Lee and Fischer [25| . 
which is a refined variant of the generic force-based model of He, Shan and Doolen [24[. The Lee-Fischer model is 
based on a square-gradient free energy functional, thus ensuring accurate reproduction of thermodynamics [60| . It 
differs from the original He-Shan-Doolen model in that it employs a special discretization scheme for the derivative 
operators and uses an alternative, but thcrmodynamically equivalent way to compute the interaction force. Since the 
effect of the discretization generally becomes noticeable only at large wavenumbers, it is expected that both models 
are equivalent in the low-wavenumber region. 

It must be remarked that the Lee-Fischer model shows a spurious drift (typically, an increase) of the total mass in 
a simulation box due to non-conservative discretization of the derivative operators (6~H . |65| . The magnitude of this 
effect is found to depend on the strength of the velocity gradients that exist in the system. If not properly handled, 
results for the structure factor can be spoiled seriously. A possible way to enforce mass conservation is to rescale all 
populations fa after each timestcp by a global factor x that is computed from the overall mass increase in the system, 
x = J2 r Am(r, t)/J2r TO ( r ' 0)- We found that this procedure gives best results for the structure factor. Despite these 
insufficiencies of the Lee-Fischer model, its behavior under thermal noise is nevertheless worthwhile to investigate, 
as this model has the advantage of not being plagued by spurious momentum currents [25l |62| and thus provides a 
promising candidate for the simulation of a range of complex fluid problems. 

In our non-ideal fluid simulations, we employ a simple Landau double-well free energy density j25l |66| 

Mp)=(3(p-Pv) 2 (p~Pl) 2 , (70) 

where pl, Pv are the prescribed equilibrium liquid and vapor densities and /3 is a compressibility parameter. The 
associated equation of state is given by po = pd p fo — fo, from which the speed of sound follows as c 2 = dpo/dp = 
pd 2 fo/dp 2 . Typically, instead of (3, one rather prefers to specify the speed of sound, which is related to (i by 

o = c 2 (p ) 

2pa(p L - pv) 2 ' 

where po denotes the reference density for which c s was computed. The shape of the equilibrium density profile 



2 Note that our variable fi corresponds to m v in ref. Il6ll 
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FIG. 1: Thermal noise in the force-based model for n = 0.08, c s — 0.265, r = 1.0. Equilibration ratios of (a,b) the density, 
(c,d) the momentum, (e) the transport and (f) the ghost modes. In (b) and (d) the dependence of the equilibration ratio of 
the density and momentum on 8, when k = (cos 9, sin 9) k, is shown for several magnitudes of k. In the remaining plots, each 
data point represents an average over all directions in fc-space at each magnitude \k\. j x denotes the rr-component, ju the 
longitudinal and j t the transversal component (with respect to k) of the momentum j. 



assumes (far from the critical point) the mean-field form 



. . 1 . . 1 . . 2y 

PMF(y) = ^{PL+Pv) + 7,KPl -pvjtanh — , 
2 2 w 



where the interface width w is given by 



Ik 1 



w 



P PL" PV Cs 

with k being the square-gradient parameter. Finally, the surface tension can be expressed as 

_ {PL-PV? /7T-3 



(71) 



(72) 



(73) 



A. Equilibration tests 

In a linear, homogenous fluid, Gaussian thermal noise described by eq. is expected to produce Gaussian 
distributed fluctuations of the LB modes with a covariance matrix given by cq. (|5T|) . We test this basic result 
by performing simulations in a uniform one-phase system of mean density p = 1.0, choosing the square-gradient 
parameter and the speed of sound as K = 0.08 and c s = 0.265. The simulation box is a two-dimensional periodic 
domain of size 128 x 128 lattice units (l.u.). The fluctuation temperature is chosen as T = 10~ 7 (setting fee = 1 hi 
l.u.), and all relaxation times are set to a value of r = 1.0. Note that the stability constraint of LB, \u\ <C a s , together 
with relation (|39p leads to an upper bound of o\pq on the fluctuation temperature T fl8l ]. Simulation results are most 
conveniently compared to theoretical expressions by computing the equilibration ratio, which is defined as the ratio 
of the equal-time correlations (\5m a (k)\ 2 ) of a LB mode divided by its expected variance G aa (k). This quantity is 
averaged over 400 simulation snapshots. The equilibration ratio is computed for each wavenumbcr k spanning the first 
quadrant of the first Brillouin cell of the reciprocal lattice. Here, a wavevector of k a = 7r corresponds to a wavelength 
of 2 l.u. The basis set T a used to obtain the modes m a from the populations fi is presented in appendix [Bj For this 
choice, the modes are given by m a = {p, 

Jx> jyi c,p ww , Pxy, Ixi Qy-i where e denotes a bulk stress mode, p ww and p X y 
are shear modes, and q x , q y and e are ghost modes. 
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Fig. [TJ shows the equilibration ratios for all LB modes in the force-based model of Lee and Fischer. We observe 
that, while the momentum is equilibrated to high accuracy (Figs. [Tb,d), the structure factor (Figs. [T^,b) shows an 
unusually large error for all but the smallest wavenumbers. The longitudinal momentum mode j\\ and the bulk stress 
mode e resemble the error found in the structure factor to a certain extent, as is expected due to the coupling between 
these modes (l8j . We have further noted a significant r-dependence of the density and momentum equilibration ratio 
at larger wavenumbers, prohibiting in fact the effective use of the Lee-Fischer model for values of r much different 
from unity. These effects are not present for the non-conserved modes, which are found to always remain well 
equilibrated up to k ~ 2.5. Although a more detailed analysis of this model is out of the scope of the present work, 
one might attribute part of the observed discrepancies to the lack of microscopic reversibility of the forcing-term (see 
the discussion in sec. IIV|) . This is also indicated by the violation of global mass conservation caused by the particular 
spatial discretization scheme employed in the model. 

B. Capillary fluctuations 

The equilibration tests of the previous section are performed in a homogenous, and therefore, fully linear system. 
However, practical applications of non-ideal fluid simulations, including, for example, phase coexistence, usually entail 
the presence of nonlincaritics. Although the Langevin noise used in the present work was derived from linear equations, 
it is nevertheless expected to be applicable to small amplitude fluctuations around inhomogencous equilibrium states 
(see the discussion in sec. IIV[) . Crucially, then the local values of the density po and the relaxation parameters A Q 
must be used in the computation of the noise variance (|53p on each lattice site. 

As a basic example for such a fluctuating non-linear system we investigate in the following capillary fluctuations of 
a liquid- vapor interface (6(| [67j ■ We mention a number of previous simulations of fluctuating interfaces using lattice 
gas automata [H, [6t| and a fluctuating ideal gas LBM [7(j ■ Note that in the latter work, the interface was modeled as 
an elastic membrane, which is different from the present approach, where the interface is represented by a smoothly 
varying order parameter. 

Capillary fluctuations are excited by the thermal noise in the bulk [7l|, [72j and can be described in terms of a local 
height function /i(rii), where r 1 1 denotes a position in the intcrfacial plane. In the present, two-dimensional situation, 
rn = x, while we denote the perpendicular coordinate by y. In the classical capillary wave theory [73ll74l]. h is defined 
by the shift between the intrinsic density profile p; n t and the instantaneous fluctuating density profile p(r), 

p(v) = p ilLt (y - h(v n )) . (74) 

In the harmonic approximation, the static spectrum of the local height fluctuations h of a flat interface is given by 

<|Mk||)| 2 H^J, (75) 

where k|| is the wavevector in the interfacial plane and a is the surface tension, eq. (|73p . For the present case of a 
Ginzburg-Landau model with bulk free energy given by eq. (|70p . far from the critical point, we can take for p- m t the 
mean-field profile (|7ip and obtain h by fitting pj nt to the density profile obtained from simulation. This procedure 
smooths out (bulk-like) density fluctuations that are always present in the interface [75T - |77| and otherwise lead to 
deviations from the capillary wave theory predictions at high wavenumbers [l8| . 

In order to test whether the static spectrum can be reproduced by our fluctuating non-ideal fluid model, simulations 
of a liquid stripe in a fully periodic, rectangular box are performed. The extension of the stripe is taken as 1024 x 200 
l.u., and the box size accordingly as 1024 x 400 l.u. Simulation parameters are pl = 1.0, pv = 0.5, j3 = 0.1, 
k = 0.08, T = 10~ 7 and r = 1.0. The capillary spectrum is obtained, after neglecting an initial roughening period 
[6^, [6^], by averaging over 2000 snapshots in a simulation running for 10 6 timesteps, which is two orders of magnitude 
larger than the largest possible relaxation time of a capillary fluctuation in the system [78||. Crucially, as we are 
working on a lattice, k 2 in eq. (|75p has to be replaced by the Fourier-transform of the proper one-dimensional discrete 
Laplacian, 2 — 2 cos fc| | . The difference between the continuum and lattice Laplacian becomes noticeable only for large 
wavenumbers (k > 1). 

Fig.[2]shows the static spectrum of the interfacial height fluctuations obtained for the fluctuating Lee-Fischer model. 
We find perfect agreement between simulation results and the theoretical capillary structure factor ([75]) for practically 
all wavenumbers. We attribute the slight deviations for wavenumbers above k ~ 1 to the breakdown of the harmonic 
approximation on which eq. (|75"]) is based. These effects will have to be investigated in future works. We remark that 
a check of the variance of the velocity fluctuations parallel to the inhomogeneity indicated that the fluid remains well 
equilibrated within an error of a few percent. We finally remark that essentially identical results can be obtained 
for capillary fluctuations in the modified-equilibrium model fl9j if definition (1741) is used for the determination of the 
intcrfacial height instead of a simple crossing criterion, as was employed in [181 ] . 
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FIG. 2: Capillary fluctuations of a planar one-dimensional interface in the fluctuating Lee-Fischer model. The equal-time 
spectrum of interfacial height fluctuations obtained from simulation [dots] is compared to the theoretical capillary structure 
factor [solid line, eq. (|75[) ]. k denotes the wavenumber in the plane of the interface. Simulation parameters: pL = 1.0, pv = 0.5, 
P = 0.1, K = 0.08, t = 1.0, the interface width is approximately 5 l.u. 



In this work, we have demonstrated how thermal noise can be incorporated into the discrete Boltzmann equation 
for the ideal and non-ideal fluid. The DBE is a precursor of the LBE, the latter being a well-established tool for 
fluid dynamical simulations. By linearizing the DBE and promoting it to a Langevin equation, fluctuations can be 
analyzed within the theory of non-equilibrium thermodynamics due to Onsager and Machlup. The covariance of the 
Gaussian noise sources in the fluctuating DBE follow from a fluctuation-dissipation relation that is expressed in terms 
of the DBE time-evolution operator and the equilibrium correlations of the distribution function. The equilibrium 
correlations are determined by invoking results of the kinetic theory of fluids. 

Technically, the analysis of the DBE is simpler than of the LBE due to the linearity of the Fourier-transformed 
advection operator and the absence of lattice corrections that occur in various places of the LBE. We have shown how 
fluctuating LB models can be constructed in a straightforward way starting from the DBE: As the noise represents just 
another source term in the DBE (in addition to a possible body force) the noise covariance of the associated fluctuating 
LBE follows immediately using well-established relations between both equations. Indeed, using the present approach 
to re-derive the noise covariance of the ideal gas model and the non-ideal fluid model of Swift et al. led to expressions 
in agreement with previous works based on the LBE [lj| [IH ■ 

As a central result of the present work, we applied the theory to the DBE of a non-ideal fluid model in which the 
thermodynamic interactions arc mediated by a body force. We have outlined general requirements that such a model 
has to fulfil in terms of the structure factor and the expression of the interaction force in order to obtain a physically 
sound and mathematically well-defined noise covariance. In particular, the noise should neither be affected by forcing 
nor advection, as both terms represent reversible contributions to the dynamics. If these requirements are met, the 
noise in the LBE of a force-based non-ideal fluid is of the same form as in the ideal gas model. 

Simulations of the fluctuating LBE version of the Lee- Fischer model indicated that satisfactory equilibration at all 
but the shortest length scales is indeed possible in this model with ideal-gas-type noise. However, we also observed 
spurious influences of the relaxation time on the quality of the results. These findings are most likely the result of 
a not fully reversible forcing term of the Lee-Fischer model - which also becomes manifest through the violation of 
mass conservation - and have to be clarified in the future. 

The further development of force-based FLBEs that overcome above mentioned deficiencies is now clearly put 
as a challenge for future works. In this regard, it would be desirable to test the present theory with other types of 
force-based models, such as two-distribution- function models derived from the original He-Shan-Doolen approach [79| . 
Furthermore, it would be extremely interesting to assess the abilities of the FLBE under non-equilibrium conditions, 
such as a fluid under uniform shear (Irl [80j. Also, the extension of the presently available "athermal" FLBEs to 
include energy conservation remains an important objective for the future. 



VII. SUMMARY 
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TABLE I: Basis set of the D2Q9 model used in the present simulations. T a i denotes the basis vector, N a its length, m a is 
the designation of the corresponding moment and A a denotes its eigenvalue in the relaxation operator. m„ q = r a i/ a eq is the 
expression for the corresponding moment of the ideal gas equilibrium distribution. 
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Appendix A: Non-ideal fluid model 

We briefly review the force-based non-ideal fluid model of He, Shan and Doolen [24[ and show that it fits the formal 
structure described in section HV B II In the model of Hc-Shan-Doolcn, the interaction force is derived from a Vlasov- 
typc mean- field approximation taking into account an Enskog volume-exclusion effect. The resulting expression for 
the interaction force can be written as [13, [zll 

F int = K pVV 2 p-V(p -p^ 2 ), (Al) 

where k is a constant ("square-gradient parameter") and po is a non-ideal pressure defining the equation of state. 
After Fourier transforming and writing 5po = c 2 s 8p 1 the linearized interaction force becomes 

6F^=ik[cl(k)-a*]6p, (A2) 

where we have used the definition of the generalized speed of sound, 

c2(k) =cl + n Po k 2 . (A3) 

In order to compute the structure factor, we note that the interaction force (|A1|) can be derived from a square-gradient 
free energy functional of the form [8l| , 

F[p\ = J dv(f Q (p) + ^\Vp\ 2 ) , (A4) 

with /o being a bulk free energy that is related to the pressure by po ~ pd p fo — /p. The structure factor obtained 
from the linearized free energy functional has the usual Ornstein-Zernike form [36ll37j. 

with c^(k) being the generalized speed of sound given by eq. (|A3[) . 



Appendix B: Basis set 

Table U shows the basis vectors T a i [l(| and the associated modes m a of the D2Q9 model used in the present 
simulations and in the theoretical analysis of the modified-cquilibrium model (section IIVB 2\ . Suitable basis sets in 
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three dimensions can be found, for example, in [ly, |42j. The parameters A a denote the eigenvalues of T a under the 
relaxation operator A of the DBE. The first three rows cover the conserved hydrodynamic moments, the next three the 
non-conserved hydrodynamic moments, and the last three the ghost (or kinetic) moments. The moment e describes 
a bulk stress mode, which is related to the deviatoric stress S a p defined in eq. (fTT)) by e = 3Tr S . The eigenvalue A e 



of e is related to the bulk viscosity Vb by v b = cr^/X e . The quantities p ww = S xx — S yy and p xy = (S xy + S yx )/2 are 
shear modes, with a common eigenvalue X s related to the shear viscosity v s by v s = a 2 s /\ s . The ghost sector finally 
consists of a ghost density mode e and ghost vector current q a , with eigenvalues A e and A g , respectively. Using the 
numerical expressions for the lattice velocities, the transformation matrix can be written in compact form as 



T = 



( 1 




-2 




V 2 



1 1 111 1 l\ 

0-101-1-11 
1 0-111-1-1 
1 1 1 4 4 4 4 
-1 1 -1 
1-11-1 
-10 1 2-2-22 
3-10 12 2-2-2 
-4 -4 -4 -4 8 8 8 8 



(Bl) 



Appendix C: Advection operator 



It is instructive to write down the explicit expression for the DBE advection operator in moment space, A a t,(k) 
-iTajk • CjT^ . Taking the basis set T a defined by eq. (|B1|) . we obtain 

/ 



A ab (k) 



k x ky 

ky 


6 2 

ky k y U 
fi 2 ^ 


. . . \ 


w ^ik^ ^iky 
2/c x 2k y 
3 3 
k y kx 
3 3 




&x ky 
kx k y 
,3 3 

ky kx 
3 3 




kx r)L 

6 2 Ah, y 

ky ky C\J 

6 2 


kx_ 

ky 

6 

Ak x Ak y . J 



\ 

This is identical to the 0(fc)-term in the expansion of the LBE-advection operator [II 

A^(k) = T a j cxp(-ik • C] )T^ = I + Aab (k) + 0(fc 2 ) . 
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